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1.0 DiTRODOCTlOH 


Th» in-traok error resulting in oMiaiting elliptic <»‘bits with Cartesian ooordi* 
mtes is reduced by several orders of magnitude with the use of a time elmmt. 

A general time element, to be used with uiy arbitrary independent variable, is 
proposed. 

Taking the exaoq>le of a transfer orbii tosr a geosynchronous mission, a eon^ari- 
sm with ttM eooentrio, te*ue, and eUiptio anomaly as the Independent variable 
is presented. The elllptio anomaly (a new anomaly intermediate betwem t)w eo- 
Mntrio and the ^ue anomaly) is ^own to perfcwm significantly better Umn the 
othM* olassioal anomalies. 
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2.0 EFPICIEHT INTEQBATION OF ECCEMTRIC ORBITS 

nils doounent stall denote, eooentrlo orbit, an elliptic orbit of 
eceentrioity e > 0.1. An «i'bit is highly eocentrio if its eooentrioity e > 0.7. 
The transftt* ca*bit linking a near-Earth, near-oiroular parking (»*bit to a near- 
olroular geosynchronous <x'blt is an exang>le of a hi|^ly eeoentrlo orbit. 

Efficient integration of an a?bit is defined as the manerioal integration of an 
(s*bit with a given accuracy and the miniatiuB mabw of steps per revolutiofi. 

In a preceding document (ref. 1) dedicated to the numerical IntegratiMi of 
highly eccentric o?bits, the following observatlms were made: 

a. An eccentric orbit should not be Integrated with equal stepsizes in time; 
l.e., time stepsizM should be small near perioenter and large near 
apocenter. 

b. An Integra tioi method with automatic stepsize control is adequate, but inef- 
ficient, tor the integration of eccentric orbits. 

0 . Better efflcimcy is obtained by using an independent variable other than 
time, by providing an analytical stepsize regulation, and by choosing a 
fixed-step integrati<xi method. 

d. For larj^ stepsizes, the iiq>ortance of a proper formulatiof: of the differen- 
tial equations of moti<xi (DEM) increases. 

e. For a certain range of stepsizes, the straightforward fonmilatlon of DB4 in 
terms of Cartesian coordinates, together with the use of a tioM element, com- 
pares favorably with more elaborate transformations of DEM. 

Automatic stepsize integration methods are based on the consideration of an esti- 
mation of the local truncation err<m* per step. Such a method does not provide 
a stepsize distribution along the orbit as regularly and efficiently as with the 
help of an analytical stepsize regulation. In addition, by having equal steps 
in tee indepmdent variable, a hi^ly efficient multistep integrati<m method can 
be used. 

Observaticm (d) is trivial for very large stepsizes. With one. step per revolu- 
titti, for instance, the Integration of Cartesian coordinates produces nonvalid 
results vdiile the integration of a set of orbital elements correctly reproduces 
the two-body part of the motion. 

For reasonably small stepsizes (> 50 steps/rev) , and for moderate propagation 
time (<50 revolutions) , a formulation of the DB4 in terms of Cartesian coordi- 
nates is adequate and gives fairly accurate results concerning the orbital 
shape. However, depending (xi the physical time, the in-track error can grow 
outside any limit. This is because of the well-known Instability of the DEM di- 
rect formulation of the two-body problem (refs. 2 and 3). 
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If Carteslixn coordinates are to be used, a special treatment of the In-track 
error Is required. This Is achieved by the introduction of a so-called time 
element . 

A time element, to be used with Cartesian coordinates, has been proposed only In 
the case of the eccentric anomaly as the independent variable (refs. 2 and 4). 

As shown In the following sections, the eccentric anomaly does not lead to the 
most efficient analytical stepsize regulation. The purpose of this document is 
to propose a general time element, which is valid for any type of anomaly. 

From a mathematical standpoint, there is little doubt that a formulation of the 
DEM in terms of orbital elements is always preferred to Cartesian coordinates. 
From a practical standpoint, however, it is often observed that the analyst 
who prepares the software for space missions is hesitant to use variables differ- 
ent. from Cartesian coordinates. This is because the observations are physical 
lengths and velocities directly expressed in Cartesian coordinates; therefore, 
a VE\luable contribution is made if the advantage of using orbital elements can 
be preserved, to a certain extent, in a Cartesian coordinates formulation. This 
is the primary concern in a time element search. 

The time element development in arbitrary independent variable terms was made 
possible by adopting the Hamiltonian mechanics framework in the extended phase 
space. Section 4.0 outlines the derivations and results. By combining the time 
transformation to be used with Cartesian coordinates (sec. 3.0) with the gen- 
eral results of section 4.0, a general time element will be derived in 
section 5.0. 

Finally, a numerical comparison in section 6.0 will show the time element bene- 
fit and the outstanding performances of the elliptic anomaly as the independent 
variable. 


3.0 TIME TRANSFORMATIONS 

An analytical stepsize regulation is achieved by a change in the independent var- 
iable (ref. 5) 


dt = f dx 


(3.1) 


where t is the physical time, x is the new independent variable (sometimes 
called fictitious time), and f is a function of the state. Equation (3.1) is 
called a time transformation. The function f is chosen to be of the form 


f = Cqi r« 


(3.2) 


where r is the radius vector. 


a is a real number, and c^^ 


is a parameter. 
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Table I indicates the four classical choices for a new independent variable, 
which corresponds to the four classical anomalies (mean, eccentric, true, and 
elliptic). The definition of the parameter involves the semlmajor axis 

a, the universal gravitational constant times the mass of the central 
body p, the eccentricity e, and the complete elliptic integral of the first 
kind K(e). The development of the formula can be found in reference 5 or, 
following a more general point of view, in reference 6. 


TABLE I.- THE FOUR CLASSICAL TIME TRANSFORMATIONS, 

f = dt/dT = Cq, 


a 


ca 


Name of angular variable 


0 

1 

2 

3/2 


Va3/p 

[uaU-e2)] -’/2 

2K(e) 

Ti Vpd+e) 


Mean ^momaly 
Eccentric anomaly 
True anomaly 

Elliptic anomaly 


The case of the elliptic anomaly (a particular choice of an intermediate 
anomaly, intermediate between the eccentric and the true anomaly) is emphasized 
here because the elliptic anomaly is, from a mathematical standpoint, the natu- 
ral way to generalize the circular polar angle to the ellipse. 


4.0 CANONICAL ELEMENTS IN TERMS OF AN ARBITRARY 
INDEPENDENT VARIABLE 


In reference 5, a set of canonical orbital elements of the Delaunay type in 
terms of an arbitrary independent variable t is proposed. The equations of mo- 
tion in the extended phase space are 



(4.1'. 


( 4.21 
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where 

bT = a, g, h) 


Is the quadrl vector of canonical elements, and 


bT = (L, y, G, H) 


the quadrivector of the conjugate elements: 

SL = time element L = total orbital energy 

i|; = anomaly 4* = conjugate of the anomaly 

g = argument of peri center G = total angular momentum 

h = ascending node H = component of angular momentum 

perpendicular to the equator 

V is a perturbative gravitational potential, and P is a trivector containing 
all perturbative forces that are not derivable from a potential. Ihe partial de- 
rivative of a scalar with respect to a vector is considered to be a column vector. 

The Hamiltonian F 


F S y - p/( 2L)^/2 + fv = 0 (4.3) 

vanishes in the phase space, where f is the time transformation (eq. (2.2)). 

The physical time t is obtained via the time element I and other canonical 
variables through the time element equation 


t 


il 


y 

+ — (u - e sin u - 
2L 


4 ') 


(4.4) 


where 

u is the eccentric 

anomaly. 




Close 

attention should now 

be given to 

the 

equations of the pair (H, 

L) of conjU' 

gate 

variables 







3f 

3r 

/3V 

\ 

3t 



V = p/(2L)3/2v — 

+ f- 1 


P) 

L' 

(4.5) 


3L 

3L 

\ 9 £ 

/ 

3L 
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where ' is an abrevlation for the derivative with respect to T. The par- 
tial derivatives of f, r, and t relative to L are 


3f 

— = cf/(2L) 


where 


c = -3f true and intermediate anomaly, 
c = -2, eccentric anomaly; 

3r 

— = -r/(2L) 

3L 

3t t 

— = - — (u - e sin u - 4») . 

3L 2l2 


5.0 A GENERAL HME ELEMENT 

The sum of two terms is recognized in equation (4.4) as 


t 


4 ' 

— (u - e sin u) 
2L 



The first term is proportional to the mean anomaly and the second term is the 
time of pericenter passage. This last term is chosen as a new dependent variable 
q , where 


q = J, _ — ij) (5.1 ) 


6 


80PM2 



The equation for t now reads 


¥ 

t = — (u - e sin u) + q (5.2) 

2L 


By differentiatiai of equatiai (5.1), the differential equation for q is 
obtained by 


q' = i' 



♦ \pL' - — 

2l2 2L^ 


(5.3) 


By assuming that all perturbations, whether conservative or not, are Included 
into the vector P, V s 0, and by using equation (4.3) the expression for 
4' is 

4'=p/(2L)1/2 

and V * 

4'* = -wL'/(2L)3/2 
Equation (5.3) now becomes 


q' = S,' ♦ 3 m\ 1/L'/(2L)5/2 _ y^'/(2L)3/2 


The equation for is needed. According to equation (4.1) and because V s o. 


9r 3t 

\|)' s 1 - f— P L' 

34'” 34' 


By using two-body relations, after some algebr? 


3r cos u 

— (2L)1/2p 

34* ue 


II 
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and 

3t 

— = (u - e sin u)/(2L) 

3Y 

throu^ (5.2). 

By taking equation (4.5) into account, q* becones 
q' s f(1 - cos u/e)PT/(2L) 

(5.5) 

♦ y/(2L)5/2 pu ♦ sin u [cos u - 2(e ♦ I/*)]? L* 

The elimination of the variable i|i, which was the purpose of defining the new 
time element q, is apparent. 

All quantities in equation (5.5) can be estimated in terms of the Cartesian 
state (r, v) with the help of the two-body relations 


r*v s ye sin u/( 2 L )^^2 


and 


r = y(l - e cos u)/(2L) 


Therefore, equation (5.5), together with equation (4.6), are ideally suited for 

supplementing the Cartesian equations of motion with a general time equation. 

A few remarks about equation (5.5) can be formulated as follows: 

a. Equation (5.5) is valid for any type of independent variable. 

b. Equation (5.5) is singular for vanishing eccentricity. This is a result 
of the development from the Delaunay el«nents, which are singular for zero 
eccentricity. However, this restriction is irrelevant because near-circular 
orbits do not require a time transformaticxi. 

c. A secular term is present in the coefficient of L' . This term say cause a 
reduction of accuracy for Icwig-term orbit propagation; however, in this case, 
using a formulaticm with Cartes! coordinates is not recommended. Instead, 
a set of elements should be used. 
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d. Vftien the independent variable is the elliptic anomaly, the only elliptical 
function to be computed tor estimating the right-hand side of the differen- 
tial equations is K(e), whi<A is a complete elliptic integral of the first 
kind appearing in the time transformation evaluation f. As discussed in 
refo'ence 6, a fast algorithm is available for this integral so that the 
corresponding overhead is of negligible effect <»i the computing time. 

To be complete, the Cartesian equations of moticm, to be integrated with equa- 
tions (4.6) and (5.5), are recalled as 


£' = fv 

v' = f(P - yr/r3) 

6.0 NUMERICAL COMPARISON 

In order to have a preliminary assessment of the efficiency of the time element 
proposed in the preceding section, a numerical comparison is made. As a sample 
problem, a standard tranaftt* crblt between a near-Earth parking orbit and a geo- 
synchronous orbit is chosen. Orbital elements of the transfer orbit are as 
fo 11 ows : 


semi major axis 

a 

r 

24 371 km 

eccentricity 

e 

= 

0.73 

inclination 

i 

= 

300 

longitude of ascending node 

n 

= 

oo 

argument of perigee 

u 

= 

2700 

initial true anomaly 

4» 

= 

0° 


As an integration method , the standard fourth-order Runge-Kutta scheme is cho- 
sen. By no means should this imply that such a simple method is adequate for 
orbit computation; a hi^er-order Runge-Kutta or a multlstep method would be far 
more efficient. However, for comparison, the primitive character of the 
fourth-order Runge-Kutta scheme will emphasize the stabilizing role of a time 
element. 

The fa"ce models, as defined in reference 7, are chosen. The following four 
cases are Investigated. 

a. No perturbations (table II) 

b. Earth oblateness perturbation (table III) 
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c. Eau'th oblateness and Moon perturbation (table IV) 

d. Earth oblateness and atmospheric drag, and Moon perturbation (table V) 

One revolution Is Integrated with four different stepslzes (25, 50, 100, and 200 
steps/revolution) and three anomalies: eccentric (oi = 1), true (a =2), and el- 
liptic (cx = 1.5), with and without time elements. 

The reference orbits were estimated by a more refined Integration. The result 
of the experiments, displayed In tables II through V, show the accuracy reached 
at the end of an Integer number of revolutions by comparlscxi with the reference 
trajectory at a given time corresponding to perigee pass. The error measured Is 
mostly In-track. 

It should first be noticed that the use of the mean anomaly for such an orbit 
would give totally meaningless results. 

With a simple time transformation without the use of a time element, tables II 
through IV show that meaningful results are obtained. An accuracy of less 
than 1 kilometer afver one revolution can be obtained with the true and elliptic 
anomaly and 200 steps/rev. 


TABLE II.- ACCURACY (KM) AFTER ONE REVOLUTION ALONG A TRANSFER ORBIT OF 
ECCENTRICITY 0.73 WITH (WITHCKIT) USE OF A TIME ELEMENT - NO PERTURBATIONS 


Step/revolution 

a 

25 

50 

100 

200 

1 

4.4 

1.2 

0.094 

0.0062 


(12 000) 

(840) 

(50) 

(3.0) 

2 

15.0 

1.2 

.079 

.0050 


C 000) 

(72) 

(4.7) 

(.30) 

1.5 

2.5 

.1 

.0058 

.00034 


(2 700) 

(170) 

(11) 

(.66) 
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TABLE III.- ACCURACY (KM) AFTER ONE REVOLUTION ALONG A TRANSFER ORBIT 
ECCENTRICITY 0.73 WITH (WITHOUT) USE OF A TIME ELEMENT - 
OBLATENESS PERTURBATION 


Step/revolutioi 



a 

25 


50 

100 

200 

1 

3.6 


1.2 

0.087 

0.0054 


(12 000) 


(840) 

(50) 

(3.0) 

2 

15 


1.2 

.079 

.0051 


(1 000) 


(720) 

(4.7) 

(.30) 

1.5 

2.3 


0.097 

.0054 

.00044 


(2 700) 


(1 700) 

(11) 

(.66) 


TABLE IV.- ACCURACY 

(KM) 

AFTER ONE REVOLUTION ALONG A TRANSFER 

ORBIT 


OF ECCBJTRICITY 

0.73 

WITH (WITHOUl) 

USE OF A TIME ELEMENT 



OBLATENESS 

AND THIRD-BODY 

PERTURBATIONS 



Step/revolution 


a 

25 

50 

100 

200 

1 

3.4 

1.2 

0.090 

0.0061 


(12 000) 

(840) 

(50) 

(3.0) 

2 

15 

1 .2 

.078 

.0045 


(1 000) 

(720) 

(4.7) 

(.30) 

1.5 

2.2 

.089 

.0045 

.00014 


(2 700) 

(1 700) 

(11) 

(.66) 
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TABLE V.- ACCURAO: (KM) AFTER 10 REVOLUTIONS ALONG A TRANSFER ORBIT 
OP ECCENTRICin 0.73 WITH (WITHCRJT) USE Cff* A TIME ELS1ENT - 
OBUTENESS, DRAG, AND THIRD-BODY PERTURBATIONS 


Step/revoluti(»i 

g 

50 

100 

200 

1 

50 

2.8 

0.17 


(15 000) 

(850) 

(41) 

2 

11 

0.74 

.047 


(250) 

(32) 

(2.5) 

1.5 

3.2 

.19 

.011 


(2 000) 

(1 200) 

(7.0) 




By including the time element, the accuracy is Improved by three to four orders 

of magnitude, as shovm in tables II through IV. For 50 steps/revolution or 

more, the elliptic anomaly gives one order of magnitude better accuracy than the 
eccentric or true anomaly. 

Without time element, the true anomaly performs slightly better than the ellip- 
tic anomaly. With the time element, the elliptic anomaly performs considerably 
better than the true anomaly, showing that the combinatlcxi time element and el- 
liptic anomaly as Independent variable reduces the In-track error most 
efficiently. 

The eccentric anomaly performs poorly without the time element, vdiile it per- 
forms as well as the true anomaly with the time element. 

The inclusion of the oblateness (table III) and third-body perturbations 

(table IV) does not seem to significantly affect the accuracy after one revolu- 
tion. This illustrates the well-kr»wn fact that the direct integration of a 
Keplerian orbit with Cartesian coordinates is unstable. Therefore, the 
perturbative effect of the force model does not add siffiificantly more 
instabilities. 


This is no longer true with the atmospheric drag perturbation. This perturbation 
acts an a hi^ly eccentric orbit like a shock at pericenter. Table V, showing 
a comparison including oblateness, drag, and third-body perturbations along 
10 revolutions, indicates a strong decrease in accuracy in all cases. The 
Integration with 25 steps/revolution had to be dropped. 

However, even in this worst-case example, a two- to three-order magnitude 
improvement in accuracy is shown by including a time element. The ellip- 
tic anomaly leaves an error after 10 revolutions 5 times smaller than the 
true anomaly, which is about 5 times smaller than the eccentric anomaly. 
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7.0 CONCLUSION 


A general time element valid with any choice of the Independent variable is 
proposed together with the use of Cartesiam coordinates for the integration of 
the elliptic motion. 

This time element is derived from a set of canonical elements of the Delaunay 
type developed in the extended phase space, and is valid for an arbitrary inde- 
pendent variable. 

By neglecting to separate the conservative from the nonconservative 
perturbations, a very simple expression for a general time element differential 
equation is found. 

Numerical comparisons show that an improvement of about three orders of magni- 
tude in accuracy can be obtained by introducing a time element for integrating 
eccentric orbits. The accuracy obtained with the elliptic anomaly (an anomaly 
intermediate between the eccentrlj and the true anomaly) can be as much 
as one order of magnitude better than the eccentric or true anomaly used as inde- 
pendent variable. 

The spectacular nature of these results, to be confirmed by more refined numeri- 
cal experiments, renders the use of the eccentric or the true anomaly obsolete 
as the independent variable when used with coordinates, while the elliptic 
anomaly combined with a time element appears to be the ideally suited independent 
variable for elliptic orbit computation. 
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